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ABSTRACT 


We have computed a stable, eq,uilibrium solution for convec- 
tion in a self-gravitating sphere of Boussinesq fluid by using a 
modal analysis in which the 6 , (p dependence of the fluid is expanded 
in the set of 168 spherical harmonics, with Z ^ 12. To com- 

pute the numerical solution of our hierarchy of nonlinearly coupled 
equations, we have developed a new relaxation method. For a Rayleigh 
number that is 30 times critical and a Prandtl number of 10, 
the flow has two orthogonal planes of reflection symmetry. The 
temperature, velocity and convective flux of the fluid as well 
as the kinetic and thermal energy spectra as functions of wave- 
length are computed. The spectra are found to be in agreement 
with both experimental observations and analytic scaling laws. 

We examine the dynamics of the energy cascade by computing the 
ratio of the amount of energy dissipation at a particular wave- 
length to the amount of energy produced at that same wavelength, 
we find that there is only a slight cascade of kinetic energy 
to smaller wavelengths but a large cascade of thermal energy. 

The fraction of the convective flux that is carried by each 
wavelength and the degree of anisotropy associated with each 
lengthscale are also determined. The stable, steady-state con- 
vecting fluid has a rotation law such that the angular momentum 
per unit mass as a function of radius is constant. 

Subject headings: cc’.vection - stars: interiors - hydrodynamics 




I. INTRODUCTION 


Thermal convection is perhaps the most important unsolved 
problem in stellar structure. In 1916 Lord Rayleigh determined 
the necessary conditions for convective instability using linear 
theory. Since then, people have been trying to calculate stable, 
equilibrium solutions to the nonlinear equations of motion. One 
of the most successful methods is finite amplitude theory (Malkus 
and Veronis, 1958), which can be used to predict the convective 
flux and temperature gradient for mildly nonlinear, low Reynolds 
number flows. Recently, Busse (1975) has computed time-independent 
convective patterns in spheres and spherical shells using this 
method. Lorenz (1963) considered the self-interaction of a single 
mode whose radial and horizontal structure was fixed for all time 
but whose amplitude was allowed to vary. He was able to compute 
a seemingly "turbulent” time dependence for the convective flux. 
Toomre et al (1977) have also computed nonlinear solutions by 
using a single mode whose horizontal structure is fixed but whose 
vertical structure is determined by the equations of motion. They 
have been able to calculate convective fluxes that have 
been applied to the second convection zone of an A star (1976). 

Each of these nonlinear methods has produced interesting results , 

but they are all, obviously, limited to low Reynolds number flows. 

A severe problem inherent in each method is that it 
requires a priori knowledge of the horizontal structure of the 
flow. It cannot be determined with these methods which, if any, 
of the infinite number of possible equilibrium solutions is 
stable. It is clearly desirable to compare the results of these 
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nonlinear analyses with the solution of a high spatial resolution 
numerical computation. 

In this paper we numerically determine the convective flow 
of a self-gravitating sphere of Boussinesq fluid for small 
Reynolds and Peclet numbers. We use the 3-dimensional modal 
analysis that was developed by this author in a previous paper 
(1979, hereafter referred to as Paper I). Although the small 
Reynolds and Peclet numbers are not appropriate for astrophysi- 
cal flows, their smallness enables us to numerically resolve 
all of the important physical lengthscales . This paper is 
necessary as a foil for the following paper in which we cal- 
culate convection for large Reynolds and Peclet numbers, but 
modal the flow at small lengthscales. The decomposition 
of the equations of motion into modes is reviewed in § 2 of 
this paper and a relaxation method that allows us to compute 
solutions to these equations is presented in §3. With this 
Qalarkin method, we are not restricted to small perturbations. 

By using a wide range of initial conditions and showing that 
they always evolve to the same steady-state, we determine 

4 

the stable equilibrium flow for a Rayleigh number of 10 and 
a Prandtl number of 10 . 

In §3 we derive different measures that can be used to describe 
the over-all properties of this flow. The. 2 and 3-dimensional 
spectra of the kinetic and thermal energies and the convective 
flux as a function of wavelengths are calculated in terms of 
modes. We also define a measure cf the anisotropy of the flow 
as a function of wavelength. In §4 we present the results 

of our numerical calculations for our stable convective solution 
these measures. Our discussion is in §3. 


in terms of 


II. MODAL EQUATIONS FOR CONVECTION 


A. Galerkin Decomposition 

The modal equations governing convection in a self-gravitat 
ing sphere of Boussinesq fluid were derived in Paper I. In 
these equations the solenoidal velocity is written as a sum of 
its poloidal part, v^, and its toroidal part, v^, which are 
derived from the scalar fields w and ip . 


V = V [3(ro))/3r] - (rV^to) e 
-p r 

v^ = rV X (\pe^) 


( 2 . 1 ) 

( 2 . 2 ) 


Each 

mean 

part 


scalar function, f(r,0,4),t) is written as a sum of its 
or horizontally averaged part , <f(r) > ,and its fluctuating 
f = f - <f> where 


<f (r,t)> 



f (r ,0 ,(j) ,t) 


/47T 


( 2 . 3 ) 


The fluctuation, f, is written as a truncated Galerkin expan- 
sion using spherical harmonics: 


(2.4) 


f ( r , 9 ,({)•, t ) - 
^cutoff 


I! I! Y'f’’‘-”C0,4,) 

2. = 1 m,Y 


wnere 


5 i 5 m _ 


,I , £ ,in _ 


2(27r)^^“ Re (Y^’™) 

m;^0 

2(7t)^^^ Y^’° 

m=0 

2(2 )^''- Im (Y^’^) 

m^O 

0 

m=0 


(2.5) 


( 2 . 6 ) 


and where Re(Y~’^) and Im(Y^’^) are the real and imaginary part 
of the spherical harmonic. The second sum in equation (2.4) 
is over OimiZ and y=R and I. The Galerkin truncation is made 
by restricting the first sum in equation (2.4) to 
We adopt the notation that <<f(r)>> is the stationary value 
of <f(r,t)>; that is , if t is a period of time that is long 
compared to the timescales over which <f(r,t)> changes, then 
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<<r 


r t 

(r)>> = j J <f(r,t)> dt 


(2. 7) 


If the value of <<f>> depends upon the origin of time used on 
the right-hand side of eq . (2 . 7), then <f(r,t)> is not a stationary 
function and <<f>> is undefined. 

We shall consider a Boussinesq fluid with a thermal expan- 
sion coefficient, a, a kinematic viscosity, v, heat capacity 
c , a thermal diffusivity, and a mean density p. The heat 


source is 


H(r) and the steady-state luminosity, i(r). 


IS 


^(r) = 4 tt J" <H 


' 2 

<H> r dr' 


( 2 . 8 ) 


VJe non-dimensionalize the equations of motion by using the 

• 3 

radius of sphere, R,as the unit of length, pR as the unit of 

mass, R'/^ac the unir of time, and £(r) /47rpc. as the unit 
of temperature. In these units the equations that govern the 
velocity, temperature, T, pressure, P, and gravitational poten- 
tial, $ , are: 


3co 


Y,£,m 


at = 


.-1 


* f’r 


- 5,(Jl+l) ^ {r e • C(vV)v]} „ 

r - - Y,Z,m 


B\p „ /dt = P (i> . ) 

, 5, ,m r Z ,Z,m 


-ZCZ+1)~^ {re . V X [vV)v]} , 

r — — Y » 5 

5 T 0 / dt = £), CT , ) 

y ,Z ,m Z y ,Z,m 


(2.9) 


( 2 . 10 ) 




JUil+1) 0<T>/9r) OJ 


Y, 


/r 


- -Cv • VT] 


Y,5.,m 


= '’r ''s 

- r““3 {r[r € • (vV)v]^, - 

r - - Ys*-»ni 


( 2 . 11 ) 


- {7 . C(v . 

= -3 


( 2 . 12 ) 

(2. 13) 


<cj> = <ijj> = 0 


(2.14) 


iill = r"^{3(r^3<T>/3r)/3r + 3^/3r 

-3 lY^r lU+1) T 

Y,S. ,m 




(2.15) 


where is the second-order differential operator defined by 
its action on f, 


(f) S [3^(rf)/3r“ - £(A+1) f/r]/r 


(2 . 16) 


The subscript, Yj stands for either i or R. The nonlinear terms 
V*C(vV)v] . are explicitly expressed in terms of u - , 

•“‘“YjXrjIU 

lb „ , and T „ in Parer I. The Prandtl nximber, Pr= v/(jL, and 

y,Z,m’ Y»^»m 

the Rayleigh number, Rs = a GR^^R) / (where G is the gravita- 

tional constant) are the two constants that appear in equations 
(2. 3)-(2.15) . 




B . Boundary and Initial Conditions 


The boundary conditions are that the surface of the convec- 
ting fluid is impermeable and stress-free. This requires that: 

= » (2.17) 

^ “ < 2 . 18 ) 

lr=l = “ 

The surface has no thermal resistivity and is therefore isother- 
mal , or 


We are also free to fix <T> to be constant for all time at the 

surface. Since the specific value of <T(1)> has no effect on 
% % 

V, T, P, oz- 3<T>/3r, we set it equal to zero. 

< T(l) > = 0 (2. 21) 


We have chosen the heat source so that it has a constant value 
for r i 0.3 and is zero outside r = 0.3. 

In dimensionless units 

^\j 

H = 0 (2.22) 

3 / 47t(. 3)^ r < 0.3 

<H>-{ Q r>0.3^ (2.23) 


This heat sou. ’'ce produces a steady state luminosity 


1 


(r) = { 


(r/ . 3 ) 


r 1 0.3 


r > 0 . 3 


} 


(2. 24) 
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The total heat flux (per unit area) at the surface is 
- 3<T> /3r and is free to vary in time. The thermal and kinetic 
energy can be stored in the fluid and r'sleased through the sur- 
face in a series of bursts rather than in a continuous flow. 
However, <<3T/3r>> = -1. 

No boundary conditions are imposed at the center of the 
sphere, but regularity requires that leading order terms of 

^,2,ra’ S.l.m 

origin. The dipole component of the velocity may be nonzero 

at r = 0 . 


III. ARTIFICIAL TIME METHOD OF SOLUTION 


A straightforward method of obtaining solutions to equa- 
tions (2.9) - (2.15) is to choose some initial data and inte- 
grate the equations forward in time. However we have developed 
a relaxation scheme, the method of artificial time, that deter- 
mines convective solutions in a more efficient manner. 

A. The Initial-Value Problem In Real Time 

Consider a fluid with its mean temperature in conductive 
equilibrium, 

r/(.3)^ r$0.3 

< T > = { . } (3.1) 

1/r r > 0.3 


and with velocity and temperature fluctuations that are ini- 
tially small. If the Rayleigh number is sufficiently large, co and 
T initially grow almost exponentially. We can understand this 
growth by computing the linearized eigenmodes to equations (2.5) 
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- (2.15); oj n and T . grow (decay) if Rs is greater (less) 

Y5JO5JI1 

than the critical Rayleigh number, ^s^^^^(£). The growth or decay 
rates are independent of y and m. The toroidal component of 
the velocity always decays initially. Rs . (Z) is given in 
table 3.1. 


TABLE 3.1 


i 

1 

2 

3 

4 

5 

6 

7 

3 

9 

10 

11 

12 1 


295 

739 

1702 

3329 

5369 

9583 

14759 

21715 

30792 

42353 

5SS09 

71.577 ! 


If Rs then the velocity decays and the mean temper- 

ature remains unchanged. (We have not found a numerical exam- 
ple of sub-critical, finite-amplitude convection). The initial 
growth rate of each eigenmode increases with increasing Rs as do 
the Jl-values of the most unstable modes. 

For example, £=6 is the most unstable eigenmode for 
8 

Rs = 10 . Although T and, w grow rapidly, <T> changes very slowly 

and its growth rate is determined by the conductive timescale, 

—5 8 

which is 'V 10” of the growth rates of T and o) for Rs = 10 . The 

horizontally averaged temperature does not begin to change appre- 

ciably until the convective flux, <Tv^> , is of the same order as 

'f 2 

the steady-state flux,oL'/r . The mean temperature gradient 
then becomes nearly isothermal (isothermal = adiabatic for a 
Soussinesq fluid). At this point, the velocity is large and 
chiefly made up of the most unstable eigenmodes (large Z) of 












the conductive temperature gradient. The nearly isothermal tem- 
Berature gradient is not in equilibrium with a velocity field 
with these characteristics . For steady-state equilibrium with a 
nearly isothermal temperature gradient, the velocity must be made 
up of the neutrally stable eigenmodes of that same temperature gra 
dient ( not the conductive gradient). The latter eigenmodes are 
characterized by smaller values of % than the unstable eigen- 
modes of the conductive gradient. Therefore the velocity begins 
cycles of decay 

and growth, with lower and lower values of I dominating the velo- 
city until a statistically steady-state is reached. The integra- 
tion of the velocity forward in time is limited by a Courant con- 
dition; a large velocity requires a small time-step. Because the 
time-step must be small and because the velocity and temperature 
must cycle through many stages until they reach a statistically 
steady state, finding solutions by integrating equations (2.S) - 
(2.15) forward in time is inefficient. 


B. The Initial- Value Problem in Artificial Time 


Instead of integrating the equations of motion in real time, 
we can force the fluid to go through a series of states in which 
the mean temperature gradient is always in instantaneous equili- 
brium with the velocity field and convective flux. Tlie only modi 
fication to the equations of motion (2.9) - (2.15) is to set 
3<T>/3t = 0 or equivalently to replace equation (2.15) with 


3<T> 

3r 

+ 


z 


Y , . , 2 

(r) / r 


Y,^,m Y,l,m ‘^y,Z,ni 


Z(Z+1) /r 


(3.2) 
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With this modified set of equations,?, u and 3<T>/3r all 
initially grow exponentially. Within a few eddy-turnover times 
the solution converges to a statistically steady state. Any 
solution to the artificial time equations (2.9) - (2.14) and 
(3.2) in which 3<T>/3r does not change in time is also an exacT 
equilibrium solution to the real time equations although it is 
not necessarily a stable solution. Stability is tested by add- 
ing a perturbation to the equilibrium solution, integrating it 
forwcrd in real time, and testing to see whether the perturba- 
tions grow or decay. Regardless of the amplitude of the pertur- 
bation, we have never found a time-independent solution generated 
from the artificial time equations that was unstable in real time. 

Sometimes the artificial time equations do not converge 
to a steady state but only to a statistically steady state (in 
which the mean temperature gradient , heat flux and spectra of 
The Jcinetic energy and thermal fluctuations are stationary). In 
these cases the variations in time of 3<T>/3r are small but they 
do not have the same time dependence that 3<T>/3r has in real 
time (see, Marcus 1980). If the statistically steady artifi- 
cial time solutions are used as the initial data in integrating 
the equations in real time, the solution converges to a new 
statistically steady state in only a few eddy-turnover times. 

IV. DESCRIPTION OF THE FLOW 

In our modal representa'^ion of the fluctuating quantities, 
we have chosen = 12, which requires 168 modes , (y , ,m) . 

By using a radial grid of 128 zones, each scalar is represented 
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by 21,5 0^'- numbers. It therefore becomes important to decide which 
quantities' should be calculated to describe the flow in a sensible 
manner . 


A. The Temperature Gradient and Mean Thermal Energy 

The rate of change of the mean thermal energy of the fluid 
is determined by integrating equation C2.15) twice over the en- 
tire radius of the sphere: 

3C / <T> dr^] /3t = 4 tt I <K> r“ dr + 47 t 1^=1 


The first term on the right-hand side of equation (4.1) is the 
rats at which energy is pumped into the fluid from the heat source; 
rhe second term is the rate at which energy is conducted away 
from the surface of the sphere. If the integration of equation 
(2.15) is left as an indefinite integral, we obtain the equation 
for the mean thermal heat flux : 


- 3 <T>/3r + r 


. -1 z 

+ r . 


Y,Jl,m '^Y,il,m ‘^Y,£,m +oC-(r)/r 


+Xc 


■f 


= -r"^ 3(/ dr'r' ^ <T>) /3t 


(4.2) 


The first term on the left-hand side of equation (4,, 2) is the con- 
ductive heat flux, and it is therefore important for us to com- 
pute the mean temperature gradient. The second term on the left- 

hand side of eauation (4.2) is the convective flux Each 

mode, (YjS-^m), contributes to the flux with no cross terms be- 


tween nodes. We will find it useful to consider the convective 

flux carried by all modes with a oarticular value of £ , F (£): 

con 


= r“^ Z(Z+1) 

con' 


with 


.-I 


•p - r { 0 ) 

^con Z "con 


I 

Y jm 


CO 




T 


(4.3) 


(4.4) 


The third term on the left-hand side of equation (4.2) is the 
heat flux of the steady- state fluid. The right-hand side of 
equation (4.2) is the rate at which the stored thermal energy 
is being released. 

B. Energy Spectra in 2 and 3 Dimensions 


The kinetic energy per unit mass in a shell of radius r 
is (see Paper I); 

(r) = Z! KS (£,r) (4.5) 


KE 


with 


KE (£,r) 




Z ,m 


r,z 


£(£+ 1 ) 




(4.6) 


where KE(£,r) is the kinetic energy due to all modes, (y ,£ ,m) 

for, all Y and m. Similiarly we can define a fluctuating thermal 
energy^ : 


(r) = 1/: 


<T'"> 


= E 


TE(£,r) 


(4.7) 


^TE(r) is usually called the temperature variance, but we reserve 
the term "variance" for another quantit}^. 
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with 

TE a,r) = 1/2 <L T _ (4.8) 

J,,m 

Traditionally one calculates the fluctuating thermal 
energy and the kinetic energy spectra as functions of a 3- 
dimensional wave numbers, whereas TE(A,r) and KE(£,r) are func- 
tions of radius and the 2-dimensional wave number, i. 

In a bounded spherical geometry it is more 

natural to compute these 2-dimensional spectra. Three-dimen- 
sional spectra are better suited to an unbounded or periodic 
Cartesian geometry. However, to compare our results with 
other published work, it is necessary to transform our 2-dimen- 
sional spectra into 3-dimensional spectra. Unfortunately, 
it is not a unique transformation. We provide an example with 
the thermal energy spectrum. 

The correlation of the temperature is defined: 


r 3 

C(r) = 1/2 / T(x) T (x + r) d X C4.9) 

The integral is over all space, but because T(x) is defined 
only for |x| - 1, it is necessary to insert a window function, 
W(x,r) , into the integral in equation (4.9). The simplest choice 
for WCx,r) is 

W(x,r) = { ^ if ix| S 1 and |x+rl i Ij ( 4 . 15 , 

0 otherwise 


There are several problems with this definition of the window 
% 

function. If T were independent of position we would expect 
the correlation function to be independent of r. Using equation 
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(4.10) we find that C(r) decreases with increasing |r|. The 
decrease is due to the finite geometry and not due to an inher- 
ent property of T. The window function must be chosen to compen- 
sate for this effect. Another consideration in the choice of 

'\j 

W(x,r) is that our integrations are numerical and T is repre- 
sented on a radial grid with a truncated expansion in 0 and 0. 
W(x,r) must be chosen so that C(r) reflects the properties 

'\j 'V; 

of T and not the numerical representation of T. 

It is convenient to have W(x,r=0) = 1 so that 

,1 

2 


C(0) = 47T J" TE(r) r^dr 
o 

We define k) to be the Fourier tranform of C(r) 

(5(k) = J" C(r) exp(-ik*r) d^r 


(4.11) 


(4.12> 


Integrating C ) over all directions, we obtain the 3-dimensional 
thermal energy spectrum that is a function of k=|kl. 


TE(k) = (2 




d() dr dx T(x)T(x+r) sin (kr) W(x,r) (4.13) 


By choosing W(x,r=0) = 1 we find that the integral of TE(k) 
over all k is equal to the total fluctuating thermal energy. 

2 


/ 


00 

TE(k)dk = 47 t / TE(r) r^ dr 

*^o 


(4.14) 


If W(x,r)=l, then TE(k) reduces to 

% % 


TE(k) = 1/2 (27 t)^J* Tj^ k^ dSl^ 


(4.15) 


where T, is the Fourier transform of T. One final caution: 
k 

since the snatial resolution is more limited in the horizontal 
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than in the radial direction, and since the horizontal res- 
olution is 2v ^cutoff should not expect TE(k) 

1/2 

to be accurate for k > ^^cutoff^^cutoff'^^^^ • Similiarly TE(k) 
will not be accurate for k < 1. 


C. 3 -Dimensional Energy Spectrum in Terms of Modes 


In terms of modes, equation (4.13) becomes 


TE(k) = (2TT)"^J‘dn„dn_ dr kr sin(kr) W(x,r) 


E I 


X r 


■1 j ^(x) - iT,. p ^(x)] 

J. K,x.,m i,x,m XX 




Z' ,m' 


^ x+r ’^x+r 




1/2 

* ^ ’■r.I'.o (|x+r|) 


(4.16) 


The spherical coordinates of the vector (x + r) are and 

We have chosen the z-axis of the (0 ,4 ) coordinate system to be 

r r 

the 0 axis. Using the identity 


/ 




= 2ir Pj [cos e^, 


(4.17) 


and assuming tnat W depends only on |x[, [r| and the angle be- 
tween X and (X + r) ,we find that the fluctuating thermal energy 
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SDectrum is 


J 


TE(k) = 2 I d(cos0_) dr dx x kr sin (kr) 


X W ( X , r , 9 , ) 

’ ’ r+x 


T 


(4.18) 


In equations (4.17) and (4.18), Pj(cos0^_^^) is a Legendre poly- 


nomial, and 


cos (0..j_„) = (x^ + |x+r|^ - r^) / (2xjr+x|) 


x+r ‘ 


x+r| = (x^ X + 2r X cos 0^)^'^^ 


(4.19) 

(4.20) 


Not suprisingly, equation (4.18) shows that TE(k) depends upon 
the product of T . (|x|)and T n ( | x+r | ) , but contains no 


cross terms between different modes. 

When the flow is isotropic the 3 -dimensional spectrum can be 


directly related to the 2-dimensional spectrum. If the correla- 

<\i 

tion C(r) between T(x) and T(x+r) depends only on jrj and if W 


depends only on lx| and lr| then equation (4.18) reduces to 


TE(k) = (2/ 


it) J " 


dr dx x“^ kr sin (kr) 


X W(x,r) 

r 




(4.21) 


Defining 




(4.22) 


dj2_Ckx) = (4/7r)l y sin (y) (1 - 1/2 y“/k x )dy 

*^o 

we find that 

TE(k) = k"^Z! r^TE(il,r) d, (kr) dr (4.23) 

I ^ 

D. Total Energy Budgets 


Multiplying the Boussinesq Navier-Stokes equation by 
the velocity and integrating over the volume we obtain the 
kinetic energy budget: 


3 Cl/2 y 

/ 


_ p 


dr^v^]/3t = 4 ttP R 

S 

- 3v. 3v. 

dr^ (^) (3-A) 

3Xj 3Xj 


I 


<v T> r dr 
r 


(4.24) 


where the first term on the right-hand side of equation (4.24) 
is KE^^ , the total rate at which kinetic energy is pumped into 
the fluid from buoyancy forces , and the second term on the 
right-hand side is KE^^^ , a negative-definite quantity that is 
the rate at which viscosity dissipates the kinetic energy. If 
the convective flux is much greater than the conductive flux 
and if the time variation of <T> is slow enough such that 

<TV > : X/r^ (4.25) 

r 


then 



>> 



r’^<T>)/3t 


(4.26) 


and 




-2 0 - 


KE. 4ttP R 
in r s 


•1 -e 


f dr= 1.892 ttP R 

j r s 


(4.27) 


In dimensionless units the rate at which mean thermal energy is 
generated is unity. From equation (4.27) it is apparent that 
in Boussinesq fluids the rate at which kinetic energy is cre- 
ated is much greater than the rate at which the mean thermal 
energy is created. This conclusion is quite different from the 
one derived from mixing-length arguments in a compressible fluid. 

In a compressible fluid these two energy input rates are nearly the 

■» 

same. It is because a Boussinesq fluid is incompressible (and 
can perform no mechanical work) that the mean thermal energy generation 
rate and kinetic energy generation rate are not equi-partitioned. 

It is often assumed that the rate at which kinetic energy 
enters the fluid is approximately equal to where 

and are the characteristic velocity and length of the larg- 

est convective eddies. Equating this expression for KE^^ with 
equation (4.27) we obtain 


V, . I (P R ) 
big r s 


1/3 


(4.28) 


This provides a scaling law for ^>j£g that can be numeri- 
cally verified. Furthermore, if most of the convective flux 
is transported by the largest convective eddy (with characteris- 
tic temperature fluctuation 'v T^^. g), then 


<T, . V^. > I X/r' 

big big 


= 1 


(4.29) 


and 






(P R 
r s 


(4.30) 


’big 


'V 


Equation (4.28) can also be verified numerically. 

By multiplying the thermal diffusion equation by T and inte- 
grating over the volume, we obtain the fluctuating thermal energy 
budget 


I 


8 [1/2 <T^> dr^=-47T <Tv > 


/ I 

<T\ 


3<T> 

3r 


r dr 


£ 


‘in. ^2 
- 47t / <VT • VT> r dr 
'o 


(4.31) 


The first term on the right-hand side of equation (4.31) is 
the rate at which fluctuating thermal energy is created from 
the interaction of the mean temperature gradient with the con- 
vective flux. The second term on the right-hand side of equa- 
tion (4.29) is TE^^^, a negative definite quantity that is the 
rate of diffusion of the thermal energy. 


E. Detailed Energy Budgets for Each Jl-Value 


Each shell of modes with the same horizontal wavelength [i.e. 
all modes (Yj^.,ni) with the same value of ll has an inflow and out 
flow of kinetic and fluctuating thermal energy. It is useful to 
determine the rata at which energy is dissipated and created in 
each particular Jl-shell. The rate of change of the kinetic energy 
in the shell is 


3 (^■TT f KE(r ,2.) r^dr)/3t = KE. (Z) - K£ ,(£) 
Jq in out 

r 1 2 

- 4tt I <v(r,2.) • [(v*V) v ]> r dr 


(4.32) 
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where KE^^(£.) is equal to the rate at which kinetic energy is 
directly fed into the 5.'*'^ shell 

KE. (Jl) = 47T P R y SUi+1) r^co 
in r s J Y i 

Y -'o'* 

and where is the rate at which the shell directly 

loses kinetic energy through viscous dissipation 




(4.33) 


KE^^^CZ) = - 4u 


iT 


“^3^CrKE(A,r)/3r^ 


(4.34) 


+ V {-Cil(Jl + l)]^ oj £). (to 0 ) 

L, y,i,m Y,^jm 


Y ,m 


-2 


-£(£+1) r”^ 3(rto . 3[r (to . )]/3r 

y 5 x» 5 m X/ y 5 X, 5 in 


r ar 


In equation (4.32), v(r,5,) is the velocity associated with the 
.th 


A““ shell at radius r. If the flow is stationary the left hand 

side of equation (4.32) is zero, but <<KE. (£)>> ^ <<KE ^ (£,)>> 
^ in out 


because -477 J <v(r , il) • £ ( v* V) v] > r"^dr, which is the rate at which 

,th 


kinetic energy cascades in or out of the Jl shell from other shells, 
is not necessarily zero. By comparing the relative values of 


KE.^(Jl), KE^^^(Jl), and 


-477^* <v(r, Jl) •C(v*V)v] > r^dr, we can deter- 
mine whether an Jl-shell is part of the production, dissipation, 
or equilibrium range of the spectrum. 

The rate of change of fluctuating thermal energy in the Jl 
shell is given by 


th 


3[4 


1477 f TE(r, 
•Jo 


Jl) r^dr]/3t = TE. (Jl) - TE ^(Jl) 

in out 


- 477^* <T(Jl,r) [(v*V) T] > r^dr 
o 


. 2 .. 


(4.35) 


Therate of fluctuating thermal energy input, TE^^(Jl),is 
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TE. (JO = 
xn 


J v.m ' 


r ^ 1’ . . ZU+1) dr 


(4.36) 


and the rate of fluctuating thermal energy dissipation is 


I 


TE^„^(J^) = -47t./- COT,, , „/3r)^+Jl(Jl+l)r"^T^„ . 1 r^dr(4.37) 

out y,m YjJ^’jHI Y,*,m. 


The rate at which fluctuating thermal energy cascades in or out 
of the shell from neighboring shells is 

-47t r^(Jl,r) C (v*V)T] r^ dr, where T(Jl,r) is the thermal fluctua- 
-*0 

tion of the Jl shell at radius r. 


F. Vorticity, Isotropy, and Angular Momentum 

The vorticity, ^ is solenoidal and can be expressed in terms 
of modes as 

= V[3(rii;)/3r] - (r'7^ip) e^ - rV [(V^oj) e^] (4.38) 


One oossible test for the isotropy of the velocity is to com- 
* * 

pute the helicity, </3.*v>, which is zero for an isotropic field. 
In terms of modes , the helicity is given by 


<vJ2> = J2-(Ji+l) r 


-2 


I 

Y,A ,m 


2Jl(Jl+l) 0) 


y,£,m 




- rijj g 3^(rw, g ^)/dr‘ 


(4.39) 


From equation (4.39) it is clear that the helicity is made up 

of terms that are proportional to ijj , co . and its derivatives . 

y 5 JO jin y 5 X/ ) jn 

We normally expect that the poloidal part of the velocity is domi^ 
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nated by the low i modes, whereas the toroidal component of the 
velocity (as well as the vorticity itself) is dominated by the 
high Z modes. Therefore the helicity is less a test of isotropy 
than a measure of the correlation between the large scale 
and small structure of the velocity. In a calculation where 
the large scale communicates with the small scale by a series 
of cascades (each of which introduces its own phase shifts) the 
helicity is not a useful measure of isotropy. Instead, we shall 
use a less .rigorous but more straightforward and revealing 
method to determine the anisotropy of the velocity. We define an 
isotropy function, I(£,r) as a function of r for each 5,-shell to be 

I0l,r) E 1/2 KEpj (Jl,r)/KEj^(£,r) (4.40) 

where KEj^(Z,r) and KEjj( 5-,r) are the contributions to the kinetic 
energy, KE ( 5, ,r) ,that are due respectively to tdie radial and horizontal 
components of the velocity. If I(Ji,r)=l, the flow is defined as iso- 
tropic. There are several advantages in defining a measure of 
isotropy as a' function of both Z and r. First of all, convection 
is not by nature isotropic; the buoyancy force is in the radial 
direction, so we expect the radial component of the 
velocity to dominate. We expect that I(5.,r) will be less than 1 
when Z is part of the production range. For larger il-shells we 
expect that the velocity field will no longer feel the direct effect 
of the radial buoyancy; as the energy cascades downward 
the velocity should lose its memory of the radial direction and 
become more isotropic. This return to isotropy can be tested 
by determining whether I(il,r) increases with increasing Z. VJe 


-25- 


also expect I(5,,r) to vary with the radius. At the outer boundary of 
the sphere, the radial component of the velocity goes to zero. 

The outer boundary-layer should be characterized by large values 
of I(£,r). We determine the beha.vior of I(£,r) by expressing it 
in terms of modes : 


I(il,r) = 


I. 


Y,m 


(il + 1) 0) 


Y,Jl,m 




(4.41) 


By using the regularity conditions at the origin, equation (4.41) 
shows that 


Lim 1(5,, r) -»■ 1 for all & 

r -*■ 0 


(4.42) 


Equation (4.42) is due to the singularity of the coordinate sys- 
tem at the origin (i.e. there is no distinction between radial and 
horizontal), so I(£,,r) is not a useful measure of isotropy at the 
origin. (However, at the origin I(5,,r) is the ratio of two 
small, numerically computed quantities. By comparing it to its 
known analytic behavior at the origin, it provides us with a use- 
ful test of the accuracy of our numerical code.) 

One more useful quantity to compute is the distribution of 
angular momentum in the fluid. All of the angular momentum re- 
sides in the 5, = 1 toroidal modes. Defining j(r) to be the 
angular momentum per unit mass at radius r we find (Paper I) : 
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- 1/2 


l(r) = 2(3) ■ r e ) (4.42) 


V. RESULTS OF NUMERICAL CALCULATIONS FOR RAYLEIGH NUMBER =10 

AND PRANDTL NUMBER =10 


A. Convergence to a Unique Stable, Equilibrium Solution 


We have integrated the artificial and real time equations for 
Rs' = 10*^ 33.8 and Pr = 10, using a code that is second- 

order accurate in both space and time. The solution to the arti- 
ficial time equations was found to be time-independent. This solu- 
tion, plus a perturbation, was used as the initial data in the real 
time equations. The perturbation decayed, and the solution recon- 
verged to the steady-state solution of the artificial time equa- 
tio,ns. We performed the integration several times: we used a radial 
grid that varied between 64 and 128 2 ones , an ^^utoff 
Galerkin truncation that varied between 6 and 12, several different 
time-steps, and a large number of different initial conditions to 
the artificial and real time equations. (All of the initial con- 
ditions had total angular momentum equal to zero.) In each case 
the integrations converged to the same time-independent solution. 
(When ^cutoff less than 12 we found that the small I modes 
did not change, but there was some change in the modes with I near 

^ 4 . jrjr- (See §5c. ) 

cutoff 

The final, stable solution has two orthogonal planes of re- 
flection symmetry. Despite our attempts to introduce pertur- 
bations without this symmetry, the solution stubbornly remained 
reflection symmetric. Modes with other symmetries, such as 
dodecahedral modes, were found to be unstable. 


A solution that is reflection symmetric with respect to both 


the x=o plane and c=o plane can be written so that co and T are 
made up only of ( R , i=even ,m=even) and (I , A=odd ,m=odd) modes and 
lb made up only of (R,£ = even,m=odd) and ( I , i,=odd ,m=even) modes. 
However, in no case did the planes of symmetry coincide with 
the x=o, y=o, or z=o planes; therefore, all modes are present 
in our final solutions. In fact, different initial conditions 
always lead to different final values of w „ , ip - , T . , 

? 5 and „ because the planes of symmetry always lie in 

different directions. We oonsider it an excellent test of our 
numerical code (especially the cumbersome nonlinear terms) that 
all of the different initial conditions lead to solutions with 
different planes of symmetry and that the solutions are identi- 
cal when they are rotated such that their planes of symmetry 
coincided with the x=o and y=o planes. 

Notice that a solution with 2 orthogonal planes of reflec- 
tion symmetry has its angular momentum identically equal to zero 
in all of its radial shells; if the solution is oriented such 
that the planes of reflection lie in the x=o planes, then 

ip-n „ -1 = > 1^0 n * n = 'I^T n T t , therefore j(r)=o for all 

^R,fi.=l,m=l ^R, £=$ ,m=0 ^I,i=l,m=l ’ 1 

r. Conservation of angular momentum requires only the integral 


:onstrai.nt / ' 


Cr)r~ dr=0. The system was not constrained to have 
j(r)"0 for all r. We also note that although the .2. = 1 component 
of the toroidal component of the velocity vanishes identically, 
the “oroidal velocity is about one-tenth the poloidal velocity 
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B. Temperature Gradient and Heat Flux 

The mean temperature gradient of the steady-state convec- 
tive solution is shown in Figure 1 (solid line) along with the 
temperature gradient of the conductive solution (dotted line). 
The convective gradient is, of course, more isothermal 
than the conductive solution, but it still shows a cusp 

at r = 0.3, which is the boundary of the internal heat source. 
Apparently the convective velocity and temperature fluctuations 
are sufficiently local that the convective flux cannot 
completely smooth out the cusp at r = 0.3. Because the convec- 
tive solution is steady state, 1/r 8<T>/3r must be equal to -1 
at the outer boundary. The slope of 1/r 3<T>/3r must be equal 
to zero at the origin. The temperature gradient is positive be- 
tween r = 0.55 and r = 0.86, which means that the convective flux 
is so large that the conductive flux is negative in this region. 
A negative conductive flux is a familiar property of convection 
in plane-parallel geometries (Herring, 1963). Figure 2 shows 

the ratio, F^^^/F™^. , where F„^, is the total flux and is 
con lotal Total 

equal tocL^(r)/r^ in steady-state convection. The ratio 

^con'^^To'*-al ratio of heat flux that is carried by con- 

vection to the total heat flux. The ratio is near unity over 
most of the radius of the sphere except in the outer boundary 
layer where it must go to zero. In §4 it was shown that the 
slope of ^QQj^/^Total zero at the origin (unless the 


r 
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leading order terms of ^ or vanish at the origin). 

The mean temperature as a function of radius (Figure 3) prominently 
shows a ’’footprint” of the r = 0.3 boundary. Boundary conditions 
necessitate that <T(1)> is zero for both the convective and con- 
ductive solution, but the value of <T> at the origin can change. 

In conductive equilibrium <T(0)> = 4.0, whereas in convective 
equilibrium <T(0)> = 0.68. 

C. Energy Spectra 

' The 2 -dimensional spectra of the kinetic and fluctuating ther- 
mal energies evaluated at r = 0.5 are shown in figures 4 and 5 . 

Both curves clearly exhibit the exponential decay that is charac- 
teristic of a dissipation spectrum. The total rate at which kine- 
tic energy is dissipated from the fluid, is determined from 

equation (4.22) and is found to be 4.96x10^. The Kolmagorov length- 
scale, n» is defined (Tennekes and Lumley, 1972) 

n = (P^^ / <<KE (5.1) 

r out 

and is equal to 0.212 in dimensionless units. 

The equivalent thermal length-scale, ri 0 > for which thermal 
diffusion becomes overwhelming is defined as 



and is equal to .067. With 128 grid points and I = 12 

CUuOff 

both of these scales should be numerically resolvable. Because 
n is greater than rig , we expect the dissipation part of the kine- 
tic energy spectrum to extend to lower values of 2 , than 
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in the fluctuating temperature spectrum. Nearly the entire 
kinetic energy spectrum is exponential wliile TE(£, r=0.5) is 
fairly flat for £ < 3. This flatness is 

due to the combined effects of production, cascade and dissipa- 
tion. Because all of the effects are squeezed into a short range 

of the spectrum, there is no clearly identif icable production or con- 
vective-viscous subrange. For £ > 3, TE(£,r) is dominated by dis- 
sipation. 

Both TE(£,r) and KE(£,r) curl upward at the high £ end of 
the spectrum. This is a typical phenomena of a Galerkin trunca- 
tion and is due to the fact that energy can cascade down the spec- 
trum until it reaches ^cutoff energy can cascade 

no further, so the velocity at ^(^^.^off forced to increase until 
^^out ^ '^cutoff ^ large enough to dissipate all of the 

accumulated energy. The fact that TE(£,r) curls up more 

than KE(£,r) is due to the fact that Pr > 1. Most of the kinetic 
energy can dissipate before cascading into However, the 

fluctuating thermal energy is not dissipated as efficiently, so more 
cascades down to ^(^^-{-off being dissipated, resulting in a 

greater curl. We have shown that the upward 

curl at is a truncation effect by repeating the cal- 
culation with £ ,^ = 10. In the latter calculation the spec- 

cutoff ^ 

tra curl at £ = 10 while with £ 4 . = 12 there was no curl 

cutoff 

at £ = 10. If the integration is performed with a Prandtl 
number of 0.1, KE(£,r) curls up more than TE(£,r). 

It is because TE(£,r) and KE(£,r) both decay exponentially 
at large £ that our modal truncation is valid. If we found that 
KE(£,r) or TE(£,r) decayed slowly as a function of £, the modal 
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truncation would not have been a good approximation. 

We have computed the kinetic energy spectrum as a function of 
the 3-dimensional wave number, k. According to Corrsin (1964) 
the 3-dimensional dissipative kinetic energy spectrum should be 
of the form 

KE(k) = exp C-3a(kri)^^^ 1/2] (5.3) 

where a is a constant on the order of unity. Figure 6 shows our 

4/3 

computed spectrum with (kri) as the independent variable. Fit- 
ting the slope of the curve in figure 6 to equation ( 5 . 3 ) between k=4 and 
k=S, we have determined a to be 3.07. While we should be some- 
what cautious in our interpretation of our deconvolved 3-dimensional 

spectrum (see i4b), it is interesting to note that our computed 

4/3 

values of log CKE(k)] are proportional to k and not k. 

Since both TE(il,r) and KE(J,,r) decay exponentially with Z, 
it is not surprising that the convective flux, v, F^^^( ), carried 
by the shell of modes, decreases with Z. Figure 7 shows the 

fraction of the convective flux carried by the Z shell at 
r = 0.5. Nearly 90% of the convective flux is carried by the 
Jl=l modes. We conclude that the amount of flux we have neglected 
by truncating at A=12 is negligible. 

D. Energy Balance 

In figure 8 we have clotted TE. (£)/TE. and KE. (il)/KE. 

in in in in 

which are, respectively, the fractions of thermal and kinetic 
energy that each &-shell contributes to the fluid. These ratios 
both decrease rapidly as Z increases x^rith 'v 93% of KE. and 
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'v 73% of produced by the £=1 modes. Because the Prandtl 


number is greater than 1, ( £=1 ) /KE^^ is greater than 


TE. (£=1)/TE. 


Furthermore, KE.„(£)/KE.„ decreases faster than 
in in in in 


TE. (£)/TE. , 
in in 


We have plotted the ratios KE , (£)/KE. (£) and TE ^(£)/ 
^ out in out 


TE. (£) in figure 9. The ratio KE^ , (£)/KE. (£) is equal to the 
in “ out in ^ 

rate at which kinetic energy is directly dissipated (not leaked 


through a cascade) from the £ shell divided by the rate at 


th 


which kinetic energy directly enters the £ shell. We would 


e^cpect this ratio to be small if £ were in a production range 


and large in the dissipative range. The fact that KE^^^C£)/ 


KEin(£) is approximately 1 for all £ can be interpreted in the 


following way: there is no (or very little) kinetic energy cas- 

th 


cade, so the amount of kinetic energy that is generated in the £ 

,th 


shell is directly dissipated by the £ shell and not transferred 
to neighboring shells. This interpretation is supported by the 
fact that the computed nonlinear terms in the model Navier-Stokes 

equation (2.7) are much smaller than the viscous, pressure and buoy- 
ancy terms. The lack of nonlinear interaction is also shown by 

the fact that the poloidal velocity (which is directly fed by 
buoyancy) is much larger than the toroidal velocity, which is fed 
energy only through the nonlinear cascade. The values of 
TE^^^( £ ) /TE^^ ( £ ) plotted in figure 9 show that there is a ther- 


mal cascade. TE , (1)/TE. (1) is less than 1 which means some 

out in 


of the fluctuating thermal energy in the £=1 shell must be cas- 


cading to higher values of £. For £ ^ 2,TE , (£)/TE. (£) is 
® ’ out in 


larger than 1 which means these modes are dissipating more energy 
than they produce. The extra energy that these shells dissi- 
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pate is provided by their neighbors . The thermal cascade is 
plausible because the nonlinear terms of the thermal diffusion 
equation (2.9) are of the same order as the dissipative terms. 

A test of the accuracy of our code is that it conserves 

energy; the comouted values of K'£ , and KE . are 4.96x10^ 

* out in 

and 4.94x10^ respectively. The computed values of and 

TE . ^ are 4.09 and 4.08. 
in 

E. Structure of the £=1 Modes 

Because the 1-1 modes dominate the flow, we illustrate the 
Jl=l radial component of the velocity, iU5.+l) r”^ “l 1 1 figure 

'X, 

10 and the temperature ,T^ ^ -j^jin figure 11. The Jl=l component 

of the toroidal velocity, ipy is equal to zero for all r. For 

these two figures, we have rotated the planes of symmetry to 

coincide with the x=o and z=o plane. Figure 10 is normalized 

by the maximum value of the A=1 radial component of the velocity, 
Jl = l 

max v^ =33.8. Notice that the velocity is a smooth function 
of radius and shows no cusp at r = 0 . 3 . On the other hand, 

Ti 1 1 does show a strong peak at r = 0.3. Figure 11 is normal- 
ized by its maximum value, max Tj i “ ^*^226. It is interest- 

Jl-1 

ing to comoare the maximum value of v with its estimated 
* r 

value of 46 derived from the scaling law (equation 4.28). The 

scaling law for the velocity is therefore 

fairly accurate for this flow. The Peclet number, in 

dimensionless units, is just equal to the velocity, so the 

Jl = i 

Reynolds number is approximately equal to max v^ 
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F. The 2-Dimensional Spectra and Boundary Layer Thickness 

2 

In figures 12 and 13 we 'have plotted KE(Jl,r)/r and 
2 

TE(^.,r)/r as functions of radius for all 1. All of the 24 

graphs have been normalized by their maximum values. Except for 
2 

S,=l, KE(£,r)/r is zero at the origin. As r increases, the curves 

remain near zero until they reach a critical value of radius 

and steeply rise to apeak. The critical value of r at which the 

curves begin to rise and the radius of the peak itself both 

increase with increasing values of 1. This behavior is con- 

sistent with the fact that KE(5,,r) “v r near the origin (Paper 

I). For Z - 2, KE(£,r)/r'^ decreases after its peak to a local 

minimum and then sharply increases at r = 1.0. This is 

because most of the kinetic energy in the first peak is due to 

2 

radial motions. As v^ goes to zero at r = 1, KE(£,r)/r decreases. 

Very close to the outer boundary, the horizontal components of 

the velocity increase rapidly in order to conserve mass flux, 

2 

and KE(£,r)/r increases accordingly. 

We define X(£) to be the distance between r = 1.0 and the 
radius of the local minimum of KE(£,r)/r and we use it as a 
measure of the boundary-layer thickness for each £. Figure 12 
shows that for £ - 6, X(£) slowly decreases with increasing £, 
and remains at a nearly constant value for £ ^ 7 . The values of 
X(£) were found to be insensitive to changes in the number of 
radial grid points or to stretching the grid near the boundary. 

We can understand the behavior of X(5,) by examining the energet- 
ics of the boundary-layer. The rate at which kinetic energy 
is dissipated from the boundary- layer associated with the £'"^ 
shell is ; 




What is the rate at which kinetic energy enters the boundary- 
layer? For large values of I, the amount of kinetic energy that 
enters the boundary-layer directly from the buoyancy forces is 
small compared to the rate at which kinetic energy is advected 
in. The rate of advection is proportional to the density of 
kinetic energy in the 5,^ -shell multiplied by the surface area 

of the boundary layer ('v 4ir) multiplied by the largest character- 
istic velocity, 

27t Cv(Z, r=l)]^ (5.5) 

Setting expressions (5.4) and (5.5) equal to each other, we see 
that v(l, r=l) cancels and that X(.i) becomes indtpendent of 
v(il,r) and i for large 1. 

The thermal energy spectra, TE( 1 ,r), display some of the same 

properties that the kinetic energy spectra have. As I increases, 

the radius at which the spectra rise and peak also increases. 

This is not surprising since both TE(Ji,r) and KE(Jl,r) behaves as 
21 

r near the origin. Unlike the kinetic energy spectra, though, 

2 2 
TE(5,,r)/r shows several peaks.. For TE(Z = l,r)/r 

(the square of the function plotted in figure 11) a large peak 
occurs at the boundary of the heat source, r = 0.3, and a smaller 
peak occurs at r = 0.88, the radius at wnich the con- 
vective flux (figure 2) drops to zero. These two peaks can be 

9 

seen in TE(Z,r)/r“ for < 6 . At = 6 the two peaks merge. 
Comparing figures 12 and 13 it is strikingly apparent that 


2 2 
TS(il,r)/r has a more small scale structure than KE(il,r)/r . At 

first this structure was thought to be numerical. Testing shows 

that it is reproducible and independent of radial grid changes. 

One reason that the kinetic energy is smoother than the thermal 

energy spectrum is that the Prandtl number is greater than 1. 

Viscosity if better at smoothing the velocity than thermal diffu- 

sivity is at smoothing the temperature. As further evidence of 

the viscosity's smoothing effect, we have seen (figure 10) that 

the £=1 component of the velocity shows no evidence of the heat 

source, while the temperature shows a marked peak. The small scale 

radial structure in T , is due to the temperature and momentum 

y 5 X- jin 

advection terms, (vVT) „ and (vVv) „ in the Boussinesq 

- Y,x,,m - - Y,J6,m ^ 

equations. Both of these nonlinear terms are a sum of many 


products of T , 


Y ' ' ,m' ’ ‘^Y 


t f n t r 


im 


, , and ilJ^ . 


" ,m’ " • 


Since 


■^Y',^',m'’ “Y’' 5 ^’',m'' '^Y ' ' ’ ' ,m' ' ' peak at different 

radii, the advection terms and T^ . ^ can have small scale radial 

y 5 X jin 


structure . 


G. Isotropy 

Our measure of isotropy, I(Jl,r), is plotted in figures 14 
and 15 for £=1 and £=2. Both curves have several grid points 
near r=0 and as r goes to zero, I(Jl,r) smoothly approaches 1 which 
shows the numerical code is accurate near the origin (see §4f ) . 

Figure 14 shows that the horizontal component of the £=1 velocity- 
goes to zero at r=0.67 and then rapidly increases in the outer 
boundary layer. The obvious conclusion to be drawn from I(Ji=l,r) 
is that the large scale structure of the convective velocity is 
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markedly anisotropic. 

We have illustrated I(5,= 2,r) because it is very similar to 
all of the I(£,r) curves with A>1. For JI>1, I(S.,r) diminishes 
with increasing radius (although it never goes to zero) and then 
increases in the boundary-layer. Because I(Jl>l,r) never goes to 
zero, we might argue that the velocity for modes with 5->l is more 
isotropic than the velocity when il=l. However, there is certainly 
no^ evidence for a "return to isotropy" in the small-scale modes . 
The lack of isotropy in the small scales is not surprising since 
there is not much of an energy cascade and the velocity 
(Reynolds number ~ 3) is far from being turbulent. 

DISCUSSION 

We have computed the convective flow with a Prandtl number 
of 10 and with a Rayleigh number that is 'v 30 times the criti- 
cal value. Despite a variety of initial 3-dimensional configura- 
tions of the flow, the fluid always evolves to the same steady- 
state with two orthogonal planes of reflection symmetry. This 
computation differs from previous studies of convection because 
the equations of motion have been used to determine the horizontal 
structure of the flow; an a priori horizontal structure was not 
imposed as a constraint to the equations of motion. We emphasize 
that to determine the stability of a convective solution, 
an examination of only the large scale modes is not adequate and 
it is necessary to include small scale modes. One large scale ^ 
equilibrium solution may be unstable with respect to another 
large scale equilibrium solution due to small wavelength pertur- 
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bations. If the expansion of the equations of motion does not 

include these small modes, both equilibrium solutions will appear 

to be stable when in fact one is unstable. In our computations, 

we have included enough modes so that the largest of the pertur- 

4 

bations has a viscous dissipation timescale that is 10 times 
smaller than the dynamic timescale. It is trerefore unlikely 
that the steady-state solution that we have computed is unstable 
with respect to any of the neglected small modes. The only other 
study of convection that has a spatial resolution equivalent to 
ours was done by Gilman (1978, and references cited therein). 

Unfortunately, it is impossible to make direct comparisons to 
Gilman’s work because his calculation is for a rapidly 
rotating shell fluid in which the coriolisis force is roughly 

the same strength as the convective buoyancy. His solutions 
favor no spatial symmetry and are time-dependent. 

Perhaps the most interesting result presented in this 
paper is that for the first time we have been able to calculate 
the kinetic and thermal energy spectra as functions of wavelength 
by solving the equations of motion alone. (For a somewhat similar 
derivation of a kinetic spectrum by use of an eddy viscosity, 
see Siggia and Patterson, 1978.) We have seen that the kinetic 
energy is so quickly dissipated by viscosity that it does not 
have time to cascade from the large modes down to smaller modes . 
Since the thermal dissipation rate is slower than the viscous 
dissipation rate, it is reassuring that our calculations show that 
the thermal energy does have time to cascade down to smaller wave- 
length before it is dissipated. 
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FIGURE CAPTIONS 

Figure 1 - The mean temperature gradients for conductive equili- 
brium (dotted line) and convective equilibrium (solid line) 
as a function of radius. The derivatives of both gradients 
are discontinuous at r = 0.3, the outer boundary 
of the heat source. The convective gradient is nearly iso- 
thermal outside of r = 0.3 and is positive for 0.55 < r 
< 0 . 86 . 

Figure 2 - The raxio of the convective heat flux to the total 
(convective and conductive) flux. For 0.55 < r < 8.6, 
the ratio is greater than 1 , which means that conduction 
carries energy downward in this region. 

Figure 3 - The mean temperature of the convecting fluid as a 

function of radius. The temperature of r = 1 is zero due 
to boundary conditions. The central temperature is free 
to vary and the convective flux has reduced it from its 
value of 4.0 in conductive equilibrium to 0.68 in convec- 
tive equilibrium. The boundary layer at r ~ 0.88 and the 
boundary of the heat source at r = 0.3 are apparent. 

Figure 4 - The kinetic energy spectrum at r = 0.5 for Reynolds 
number I 3 as a function of 1. The spectra shows a nearly 
exponential decay with no trace of a production or equili- 
brium range. The slight curl upward at £=12 is due to trun- 
cation . 
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Figure 5 - The thermal energy spectrum at r = 0 . 5 as a function 
of i for a Peclet number of 1 34 . The thermal spectrum 
falls off less steeply than the corresponding kinetic 
energy spectrum. The flatness of the spectrum at £ - 3 

is due to a combination of production and cascade of ther- 
mal energy. The upward curl at £=12 due to truncation is 
more severe in this figure than it is in figure 4 for the 
kinetic energy. 


u / 3 

Figure 6 - The kinetic energy spectrum as a function of Ckri^ 

where k is the 3 -dimensional wave number and n is the 

Kolmogorov length. By fitting the slope of this calculated 

curve for 4.0 - r - 8.0 to the theoretical dissipation spec- 

47 3 

trum KE(k) = exp[-3a(kn) 72], we have computed the dimen- 
sionless constant a. to be 3.07. 

Figure 7 - The fraction of the convective flux at r = 0.5 carried 
by all modes with horizontal wave number £ is plotted as a 
function of £. Nearly 90% of the convection flux is carried 
by the £ = 1 mode. Like the thermal and kinetic energy spec- 
tra, the convective flux decreases exponentially with 
increasing £. 

Figure 8 - The fraction of the total amount of kinetic (thermal) 
energy that is produced in the £'’^ shell of modes is shown 
by the solid (open) circles. While 93% of the kineric 
energy is produced by the £ = 1 mode, only 73% of the ther- 
mal energy is produced by that mode. The figure shows that 
the production of thermal energy is more important than the 
production of kinetic energy in high £-shells. 

. ... - * 
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Figure 9 - The ratio of the kinetic (thermal) energy that is vis- 
cously dissipated by the £ shell of modes divided by the 
amount of kinetic (thermal) energy that is produced in the 
I -shell is plotted as solid (open) circles. The ratio 
for the kinetic energy is near unity for all £• and is indica- 
tive of no kinetic energy cascade; all of the kinetic energy 
produced in the shell is dissipated there before it has 
a chance to cascade to a different £ shell. For £ = 1 the 
ratio for the thermal energy shows that more energy is pro- 
duced there than dissipated; this indicates that thermal 
energy escapes from the £ = 1 shell by cascading into higher 
£ shells. For £ > 2 the thermal ratio is greater than 1 
and increases exponentially indicating that the £'*'^ shell 
dissipates more thermal energy than it produces. The extra 
energy that is dissipated is provided by the energy cascade. 

Figure 10 - The radial component of the £ = 1 velocity as a func- 
tion of radius. The velocity has been normalized by its 
maximum value of 33.8. The velocity shows no trace 

of the boundary of the heat source at r = 0.3. 


Figure 11 - The £ = 1 component of the temperature as a function 

of radius. The temperature is normalized by its maximum value 
of .0226. The temperature has an obvious peak at the boundary 
of the heat source. 
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. 2 

Figure 12 - The kinetic energy spectra KE(£.,r)/r as functions 

of radius for £=1,12. As £ increases the spectra peak at 
increasing values of the radius. The edge of the outer boun- 
dary layer is approximately equal to the radius of the 

2 

local minim\im of KE(£,r)/r near r = 1.0. For 1 < £ i 6 the 
thickness of the layer decreases slowly with increasing £. 

For £ > 6 the thickness is nearly constant. 

. 2 

Figure 13 - The thermal energy spectra TE(£,r)/r as a function 

of radius. As £ increases, the radii of the maximum values 
of these spectra increase. For £ < 6 there is 
a secondary peak near the edge of the boundary layer. At 
£ = 6 the two peaks merge and for higher £ only one peak 
can be seen. There is more radial structure in the thermal 
than in the kinetic energy spectra because the thermal fluc- 
tuation are less easily dissipated when the Prandtl number 
is greater than 1. 


Figure 14 - The ratio of the horizontal to the radial component 
of the kinetic energy for £ = 1. The ratio is a measure of 
the large scale anisotropy of the flow. 


Figure 15 - Same as figure 15 except that this figure is for the £ = 2 
component of the velocity. The anisotropy for £ = 2 is quali- 
tatively the same for £ > 2 and there is no apparent isotro- 
pization of the velocity at small wavelengths. 
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